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We investigate the mechanism of heat conduction in ordered and disordered harmonic onedimen- 
sional chains within the quantum mechanical Langevin method. In the case of the disordered chains 
we find indications for normal heat conduction which means that there is a finite temperature gradi- 
ent but we cannot clearly decide whether the heat resistance increases linearly with the chain length. 
Furthermore, we observe characteristic quantum mechanical features like Bose-Einstein statistics of 
the occupation numbers of the normal modes, freezing of the heat conductivity and the influence of 
the entanglement within the chain on the current. For the ordered chain we recover some classical 
results like a vanishing temperature gradient and a heat flux independent of the length of the chain. 
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I. INTRODUCTION 

Fourier's law of heat conduction states that the heat 
flux jg through a medium is proportional to the temper- 
ature gradient: jq — —k'VT. In the stationary ID case, 
where the heat flux is constant, this implies a constant 
temperature gradient, at least if the thermal conductivity 
K is not temperature dependent. 

The rigorous deduction of Fourier's law from classical 
or quantum mechanical statistical physics has been an 
unsolved problem for decades. In 1955 Fermi, Pasta and 
Ulam 1] performed the first numerical calculations on 
the dynamics of harmonic chains, the simplest imaginable 
model for lattice vibrations in insulators. It was found 
that harmonic chains do not find their way to equilibrium 
because the normal modes are decoupled in harmonic 
systems 011] • In one dimensional ordered chains the heat 
flux is independent of the length of the chain and the 
temperature gradient vanishes. 

In disordered chains the conductivity is reduced due to 
the Anderson localization of most of the normal modes 
leading to a finite temperature gradient. The over- 
all resistance however does not increase linearly with the 
chain length, but is proportional to its quare root. That 
means that the specific conductivity still diverges with 
length. 

Non-integrability is necessary for the observation of 
equilibration of energy and diffusive heat conduction as 
described by Fourier's law. The problem of heat conduc- 
tivity in all kinds of model systems is still an active field 
of research. Exemplarily we mention [1, llfl] dealing with 
nonlinearities and [T]| dealing with momentum conser- 
vation. Nevertheless no model Hamiltonian system, for 
which Fourier's law could be proven rigorously, has been 
found yet. 

The one dimensional case is easiest to handle and can 
be partly justified as a model of the homogeneous three 
dimensional case. However heat diffusion perpendicular 
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to the direction of the heat flux is neglected. This can be 
compensated with self-consistent heat baths, which add 
noise and damping with zero average energy flux to every 
site of the chain. Bonetto, Lebowitz, and Lukkarinen |12l | 
find normal heat conductivity in such a model, Barros, 
Lemos, and Pereira [l^ study the manipulation of the 
heat flux in such a system by changing the masses and/or 
the onsite potentials in a chain with self-consistent heat 
baths. 



The other side of the problem is the quantum mechan- 
ical description of temperature and heat which include 
thermal fiuctuations and dissipation. To achieve this 
we follow the ansatz by UUersma UM , for a review see 
Nieuwenhuizen and AUahverdyan [l5|. The considered 
system is extended by introducing a heat bath consisting 
of many environment degrees of freedom. These degrees 
of freedom will be traced out and a Langevin equation 
is obtained. The energy transfer from the system into 
the bath degrees of freedom appears as a damping term. 
Reversely the initial conditions of the bath degrees of 
freedom appear as a noise term. 

The aim of this work is to use the quantum mechani- 
cal Langevin ansatz for heat conduction in a chain. For 
technical reasons we restrict ourselves to one dimensional 
harmonic systems without self-consistent heat baths in 
the ordered as well as in the disordered case. This ansatz 
was used by Ziircher and Talkner [13, and recently by 
Dhar and Roy [l^. Dhar and Roy apply their method 
to the self-consistent heat bath model [13] and recover 
their results in the classical limit. 



Alternative systems for the investigation of quantum 
mechanical heat conduction are systems of coupled spins 
(20I [2H . The main difference compared to harmonic os- 
cillators is the finite dimensional Hilbert space with e. g. 
only two energy levels per site for spin-i. Depending 
on the types of coupling and of the choice of parameters 
normal heat conduction is found or not. 
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II. MODELS FOR HEAT CONDUCTION 

A. Disordered harmonic chain coupled to two heat 

baths 

We consider a chain consisting of I harmonic oscilla- 
tors. Xj and Pj denote the coordinate and the momen- 
tum of the j-th oscillator. The oscillators have common 
mass M but each oscillator has its own onsite frequency 
ujj. Nearest neighbors Xj and Xj^i are coupled via the 
coupling constant fj . 




(1) 

There are two heat baths denoted by a and b which are 
coupled to the first and to the last oscillator via coupling 
constants Ck- 
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The Hamiltonian of the complete system is then 
H = Hch + Ha + Hb . 



(2) 
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Initially the bath degrees of freedom are independently 
occupied according to the bath temperatures Ta and Th. 
The normal coordinates of the isolated chain are occupied 
according to the temperature Teh- Then the couplings Ci 
are switched on at time t = 0. 



1. The equations of motion 

Regarding Xi{t) and Xi{t) as known inhomogeneities, 
the equations of motion for the bath degrees of freedom 
are solved and substituted into the equations of mo- 
tion for the chain. One gets the quantum mechanical 
Langevin equations 

1 
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where we use the Einstein notation, i. e. indices appearing 
twice are implicitly summed over. The matrix C is the 
coupling matrix of the isolated chain 
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with /o and /; set to zero. rji{t) and 7(i) are the noise 
function and the damping function and are discussed be- 
low. The matrix z/^ = 5ij {Sn + Su) connects the damp- 
ing term to the first and the last oscillator of the chain. 



a. The damping function The damping kernel has 
the form 

JV 2 

j{t-t')^y^^^cos{cuk{t-t')). (6) 

If we choose the bath frequencies tOk = fcA, with the 
level spacing A, and the coupling constants Ck according 
to the Drude-UUersma-^spectrum 15] 
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(7) 



and perform the limit N oo and A — > 0, we get the 
convenient result 7(t) = 7re^'"l*l with the Laplace trans- 
form 7(s) = j?^^. 

b. The noise functions rja and rji, act on the first, 
respectively on the last oscillator of the chain 



mit) ria{t)5ii +rib{t)6ti. 



(8) 



The noise is determined by the initial conditions of the 
respective heat bath 



TV 
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a = a,b. In contrast to the damping, the noise functions 
depend on the temperature of the respective bath. The 
random character of the noise function comes from the 
unknown initial conditions of the bath degrees of free- 
dom. As the initial conditions Xi{0) and Pi(0) have zero 
average, the average of rjj{t) is also zero. Later we will 
need the symmetrical autocorrelation function of ria(t), 
which is calculated using the initial conditions: 



Kait-t') := ^{{VcitM')) + {Vait'Mt))) 
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There is no analytical expression for this integral. It 
is nonzero even for = 0, indicating that quan- 
tum fiuctuation are never absent. In the classical limit 
coth (y 2kBT ) reduces to 2kBTa/ {hio), and in the Marko- 

vian limit F oo, K^it — t') becomes (5-like, i. e. white 
noise. A detailed discussion of quantum noise can be 
found e. g. in p^ . 

It is worth mentioning that the Heisenberg equations 
of motion for the operators are identical with the classi- 
cal equations of motion. All quantum mechanical effects 
follow from the initial conditions of the chain and the 
bath degrees of freedom in the noise function. 
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2. The solution of the Langevin equations 



with the interaction matrix 



Thanks to the linearity of the problem the equations of 
motion ([5]) can be easily solved in Laplace space, where 
the convolution turns into a product and the differential 
equation becomes algebraic. By using the La pla ce trans- 
form instead of the Fourier transform used in [l9| , we will 
be able to solve equation ([5]) for any t and not only for 
the stationary case. In Laplace space the equation reads 



pm 
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C'ijXj 



M 



^^j'-^X,{s). (11) 



We perform a coordinate transformation to the eigen- 
functions Yi of the coupling matrix C, i. e. to the nor- 
mal coordinates of the isolated chain, which are standing 
waves in the case of the ordered chain. The transforma- 
tion matrix is denoted by G: Xi — GijYj, {GCG^)^- = 
fliSij. The equations of motion then read 



Qi{0) , Guiiais) + Gufibis) 



M 
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(12) 



(s^ + r!|)5,, + ^ (Gi,Gi, + GuGij) 



(13) 



Both damping and noise act on the normal coordinates 
via Gii and Gu, i. e. their deflections at the ends, where 
the chain is coupled to the baths. 

We need the inverse of the interaction matrix B{s) for 
solving p2)) for Yj{s). The entries of B{s) are rational 
functions of s. We extract the common divisor of all 
matrix entries and end up with a matrix containing only 
polynomial entries, which can be inverted using Cramer's 
rule. The entries of ^4 = B~^ are rational functions of 
s with a common denominator. The poles lie in the 
left half of the complex plane. After performing a par- 
tial fraction expansion we transform back to time space, 
ending up with a sum of decaying exponentials. In Ap- 
pendix El we will have a closer look at the inversion of 
the matrix B(s) in the case of symmetric chains. 

Inverting equation (|12p and transforming to time space 
thus yields 



fc=i 
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(14) 



r 



with the response functions F°-{t) — GikAjk{t) for 
the noise r]a{t) and F^{t) = GikAjk{t) for r]b{t). The 
equation for the momenta Qj{t) = MYj{t) is obtained 
by differentiating. The response functions are sums of 
decaying exponential functions, e. g. 



•a Xkt 



c.c. , Re(Afe) < , a — a,b. 



(15) 

That means, any contribution from the initial conditions 
1^(0) and Q{0) vanishes with time. As mentioned above. 



the averages {r]a{t)) and {r]b{t)) are zero, so {Yj{t)) is 
zero, as well. Two-point-correlations of coordinates and 
momenta are the objects of interest. 



S. The evaluation of time dependent correlations 

With the response functions Ajk{t) and F!j'{t), and the 
symmetrical autocorrelation functions Ka(t — t') of the 
noise provided, one can evaluate any symmetrical corre- 
lation like i ({y,(t), Y, (i)}) or i {{Y,{t),Qj{t)}). E. g. 
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(16) 
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The generalization to time-shifted correlations is cum- 
bersome, but straight forward. 

In the integral there is a summation over the exponen- 
tials exp{Xk{t — t')) and exp(Afc' (t— i")) from and the 
w-integration in K{t' — t"), see ([TU|. The double time in- 
tegration over exponentials and cosines can be performed 
analytically. In the limit t ^ cxj it yields the result 

AfeAfe' -t- 



Then the w-intergations are evaluated for each pair of 
roots (A/c, Xk') and summed up. 

Beside calculating the roots A^ and the coefficients for 
the response functions F",^ the evaluation of the w- 
integrals is the main numerical work. For both tasks we 
have employed standard routines from the Mathematica 
environment 



B. Symmetric disordered chains 

In this section we consider a specialization of the sys- 
tem above, namely a disordered chain with left-right 
symmetry, which means that the Hamiltonian is invari- 
ant under the exchange X„ X/_|_i_„. This implies 
fn — fi-n and LUn — LUi+i^n- The normal coordinates of 
Hch are either symmetric or antisymmetric with respect 
to commuting left and right. In particular the transfor- 
mation matrix G and the noise response functions obey 
the following relations 



Yi even 



Yj odd 
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(17) 

Even modes thus respond to the effective noise 7]^ = 
Va + Vb and odd ones to rjo — rja — rjb- The nondiago- 
nal terms in the interaction matrix B{s), see (jl3p . are 
proportional to {GuGij + GhGij) and vanish, if the cor- 
responding modes have different symmetry. Thus Bij (s) 
and its inverse A are block matrices which do not mix 
even and odd modes, allowing to invert B separately for 
each symmetry family with separate sets of roots. 

For convenient notation we define the symmetry func- 
tion p{j) and the symmetry family cr(j) 

I p for even 
P{j) ■■= < ' ' ^(j) := I P(«) = P(j)} • 



o for K; odd 



(18) 



With Fj{s) Fj-{s) one can express equation p4p in 
the following compact form 



k£a(j) 



A,kit)YkiO) + ^A,kit)Qk{0) 



(19) 



The noise response functions have the same form as 
in (jlSp . but only exponentials with Xk from the same 
symmetry as Yj occur. 



C. Ordered Chains 

A further specialization of the disordered chain is to 
set all couplings /„ equal to a constant / and all ujn to 
ujQ. In this case the normal coordinates of the isolated 
chain are simply standing waves with anti-nodes at the 
ends. 



III. RESULTS 

In this chapter we present and discuss some numeri- 
cal results obtained with the techniques presented in the 
previous chapter. 

At first we will investigate chains without disorder, i. e. 
/,; = / and LOi — ujQ. In Sec. IIII Bl we will introduce 
disorder in the couplings fi. The onsitc frequencies uji 
are kept ordered. They cannot be set to zero because the 
translation of the center-of-mass must be suppressed. 
We set all uji to a constant ujq, meaning that our chain 
is fixed on a substrate and cannot move macroscopically. 
Another possibility, often employed in the literature [4|, 
[TtI . [TsI , is to fix only the first and the last oscillator with 
an onsite potential, corresponding to a free wire spanned 
between the heat baths. At some points we will refer to 
this model as well. 

In the numerics we work with dimcnsionless quanti- 
ties. The mass M and the onsite frequency ujq fix to- 
gether with the Planck constant H and the Boltzmann 
constant ks the units of all quantities. Thus, in the 
results frequencies are given in units of loq, energies in 
[E] = hujo, temperatures in [T] = hujo/ks, currents in 
[J] = hujf^, conductivities in [G] = ujokB, coupling con- 
stants in [/] = A/cJq, momenta in [P] = \/MhujQ and 
lengths in [X] = ^ hj (mujo). 

Furthermore we fix the cutoff F = IOwqj so the cou- 
plings / and 7 and the temperatures are the free param- 
eters. Unless specified otherwise we choose a chain with 
length / = 20. 

Except for Sec. IIII Dl we will focus on the correlations 
within the coordinates and momenta of the chain in the 
stationary regime, i. e. the initial conditions of the chain 
are irrelevant. 



A. Ordered Chains 

In contrast to chapter |TT1 where we started with the 
general disordered case and specialized the system untill 
the ordered chain, we will begin with the simplest case, 
the ordered chain, here. 
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1. Energy distribution in the normal coordinates 

As the calculation of the variances is performed in 
normal coordinates, we inspect the energy distribution 
in normal coordinates first. The energy in the nor- 
mal coordinates of the unperturbed chain reads Ei = 
\Amj{Y^) + ^(Q^)- As the coupling to the heat 
baths is rather strong, we take the coupling energy into 
account by determining effective frequencies. Therefore 
we perform a second calculation with both bath temper- 
atures set to zero, i. e. the entire system will be in the 
ground state. We then identify the groundstate energy 
with the zero point energy of an oscillator with the effec- 
tive frequency fli. Furthermore we assume that the virial 
theorem is approximately valid, although the normal co- 
ordinates of the unperturbed chain are coupled among 
each other by the interaction with the baths. Therefore 
one can calculate effective frequencies fli, using only the 
kinetic energies (Qf ) /{2M): 



{Eo)i « 2(i?o'")j — 



M 2 



Then we use the effective frequencies fli to calculate the 
occupation numbers from the kinetic energies at finite 
temperatures: 



1 



The low frequency modes with large amplitudes at the 
ends of the chain are affected the most by the coupling 
to the baths. 




(a) 



(b) 



FIG. 1: 



Occupation numbers of the normal modes calcu- 
lated from the kinetic energies. Parameters: / = 1, 
7 = 2. (a) high temperatures Ta = 5, Tt = 2, 
Iflt ~ 3.54513 (b) low temperatures Ta = 0.5, 
Tt = 0.2, Tflt = 0.450906. The center-of-mass mo- 
tion shows a negative deviation. 



In |Fig. 1| these occupation numbers are plotted. Al- 
though there is a temperature difference and a heat flux 
between left and right (see Sec. IIII A3p . the occupation 

numbers essentially agree with the Bose-Einstein distri- 

-1 
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in both cases. In the high 

temperature case (a). Tat agrees with the average tem- 
perature {Ta + Tii)/2. In the low temperature case (b). 



the temperature Tfit is closer to the higher bath temper- 
ature, giving a hint to the fact that the heat conductivity 
increases with temperature (Sec. lIII A3|) . There are some 
deviations from the Bose-Einstein distribution resulting 
from the approximations made in the calculation of the 
frequencies. In the limit of weak coupling they vanish. 

Different definitions of effective frequencies are possi- 
ble. For example one could use the potential energies 
instead of the kinetic energies and would obtain different 
results. The reason for this is that the normal coordi- 
nates of the unperturbed system are coupled via the heat 
baths. We are not dealing with independent oscillators 
in the ground state, which is reconfirmed by the fact that 

we find {Y^),{Qdo >{!)"■ 



2. Temperature profiles 

We transform the correlations of the normal coordi- 
nates back to real space and calculate the energy per 
site, splitting each spring energy to the neighboring sites 



En — 



~ fn-1 {XnXn-l) — fn {^nXn+l) 



(20) 



Numerical results are shown in |Fig. 2[ The energies of 
the first and the last oscillator are close to the thermal 
energies of the respective heat bath, and like in the clas- 
sical investigations (e. g. [3]), the temperature gradient 
vanishes inside the chain. With different coupling pa- 
rameters / and 7 one can change the behavior only very 
close to the boundaries. In the low temperature case the 
energies per site are dominated by the zero-point ener- 
gies. The energies of the boundary oscillators are ele- 
vated because their effective frequencies are increased by 
the coupling to the heat baths. 



(a) 



(b) 



FIG. 2: 



Energy per site (dots) . The bars represent the ther- 
mal energies of the respective heat baths. Parame- 
ters: / = 1, 7 = 2. (a) hight temperatures Ta = 5, 
Tb = 2 (b) low temperatures Ta = 0.5, U = 0.2 



We want to eliminate the zero-point energies and con- 
struct a temperature for each lattice site. Like in the 
previous section we determine effective frequencies with 
a ground state calculation, using Eq — Hoj/2. Then we 

assign a temperature Tr using E = ^hoj coth ^ 2k]^TR ) ■ 



FIG. 3: The temperatures Tr reconstructed from the ener- 
gies per site and the zero-point energies per site. 
The bars represent the temperatures of the respec- 
tive heat baths. Parameters: / = 1, 7 = 2. (a) high 
temperatures Ta — 5, Tb — 2 (b) low temperatures 
Ta = 0.5, Ti = 0.2 



2 

1 

Y 

. 4: The heat fiux J as a function of the couphng to 
the heat baths 7 and the coupling within the chain 
/. The dots mark the maxima with respect to 7. 
Ta^5,n = 2 



In |Fig. 3| the reconstructed temperature is shown. In 
the high temperature case (a) the result looks rather the 
same as in |Fig. 2] because the zero-point energies are not 
very relevant at high temperatures. In the low tempera- 
ture case (b) however, the zero-point energies have been 
compensated and one can see a temperature profile which 
lies between the two bath temperatures, except for the 
boundary oscillators. Again the temperature gradient is 
only an exponentially small boundary effect. Note that 
the interior temperature is again closer to the tempera- 
ture of the warm heat bath, indicating a higher thermal 
conductivity at high temperatures. 



of 7 where J is maximal, increases linearly with /. In our 
model 7max starts approximately linearly with / but falls 
behind for larger /. 

b. Heat flux and chain length In the case of normal 
heat conduction the heat flux is expected to decrease re- 
ciprocal with the chain length at fixed temperature dif- 
ference Ta — Th. In agreement with the vanishing tem- 
perature gradient inside the chain fSec. IIII A 2[) . we find 
that the heat flux does not decrease with the chain length 
for I > 5, |Fig. 5| The total heat conductivity rather than 
the specific conductivity is a constant in the ordered har- 
monic chain. 



The heat flux 



We set up an equation of energy continuity 



(21) 



by differentiating ((20)) . In the stationary case we find for 
the energy fluxes from site n to n + 1 



(22) 



This result agrees with the classical formula 'power = 
force X velocity' using 'force = fn{Xn+i — A'„)' and 

(X„P„) « A i^xl) = 0. 

Due to energy conservation must be indepen- 

dent of n, therefore we write J := J^^^^. The heat flux 
J is one scalar quantity and therefore easier to analyze 
than the temperature gradient. 

a. Heat flux as a function of the coupling constants 
We calculate the heat flux for different coupling constants 
/ and 7 and plot the heat flux over the /-7-plane, see 
|Fig. 4| In general the heat flux grows with / and 7, 
but / and 7 must match each other: For a given / the 
heat flux increases linearly with 7, passes a maximum at 
7max and then vanishes like 7~^. This agrees with the 
behavior found by Rieder et al. 0] in a classical model 
without onsite potentials. In their results 7max, the value 
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FIG. 5: The heat flux as a function of the chain length /. 

Parameters: 7 = 2, / = 1. (a) high temperatures 
= 5, Tb = 2 (b) low temperatures Ta = 0.5, 
Tt = 0.2 



c. Thermal conductivity as a function of tempera- 
ture We are especially interested in the low temperature 
regime. Therfore we vary the mean temperature with a 
flxed relative temperature difference e = {Ta — Tij)/{Ta + 
Tb) and calculate the conductivity Gth ■= J/ (Ta — Tb) for 
each temperature. 

Results are shown in [Fig. 6[ In the high temperature 
regime the conductivity is a constant like in the classical 
case. In the low temperature regime it breaks down and 
behaves similarly to the Bose-Einstein occupation num- 
bers of the lowest normal frequencies (gray line). In this 
case the degrees of freedom of the chain are simply frozen 
out. 

This behavior is typical for optical phonons. In our 
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FIG. 6: Black: The thermal conductivity as a function of 



temperature. Gray: c [exp (ftojo/fcsT) — 1] 
rameters: e — 0.01, / = 1, 7 = 2 



Pa- 



of the participation number p: 



(23) 



The center-of-mass mode, which is extended over the 
whole chain yields ^ = a fully localized mode yields 
^ = 1. Results for an unsymmetrical disordered chain are 
shown in |Fig. 7| The center-of-mass mode has — I, all 
the other modes have smaller localization lengths, which 
decrease with frequency. 



model the low normal frequencies start with the value 
Wo, which means all normal modes are frozen if fc^T < 
fiiOQ. If we set ujq to zero in the inner chain, the normal 
frequencies start at zero and we find the temperature 
dependence of acoustic phonons Gth oc T^, which was 
also found in [l^. 



B. Disordered Chains 

It is well known from classical works that disorder can 
lead to a finite temperature gradient inside the chain, 
6. g. 0, 0, @- There are different possibihties to bring 
disorder into play: 

• A common choice in the literature is to choose 
the masses of the oscillators randomly. This cor- 
responds to isotopical disorder in nature. 

• Random frequencies of the onsite potentials uji 

• Random coupling constants fi between the oscilla- 
tors of the chain. 

To simplify matters we confine ourselves to disorder in 
the couplings fi in this article. The fi are chosen from 
a Gaussian distribution with mean / and width cj, and 
a cutoff that guarantees that the fi are always positive. 
In the following we use the notation / ~ f ± <7f. Tem- 
perature profile, heat flux etc. are calculated for many 
realizations of disorder and averaged over. In the follow- 
ing (•) denotes the ensemble average. 

As indicated in Sec. Ill Bi the symmetry of the chain is 
relevant. Therefore we always compare the results of the 
symmetrical and the unsymmetrical disordered chain in 
this subsection. 



1. Normal modes in disordered chains 



20 
15 




1 1.25 1.5 1.75 2 2.25 2.5 

FIG. 7: The localization length ^ in dependence of normal 
frequencies Q for an unsymmetrical chain with fi = 
1 ± 0.2 averaged over 50 realizations, the error bars 
show the standard deviation. 

In the symmetric chain the normal modes are con- 
strained to be symmetric or antisymmetric, and are 
therefore less localized. 



2. Occupation numbers 

In |Fig. 8| the occupation numbers of disordered chains 
are shown. In the symmetrical chain the frequencies are 
smeared by disorder, compared to the ordered case in 
|Fig. 1| (a) . At the same time the occupation numbers still 
agree with the Bose-Einstein distribution at the mean 
temperature. 

In the unsymmetrical case, |Fig. 8| (b), this changes: 
The distribution broadens and the points lie between the 
Bose-Einstein distributions corresponding to Ta and Ti,. 
It is striking that the points belonging to the high fre- 
quency modes tend to lie close to either bath tempera- 
ture. This is due to the fact that the strongly localized 
high-frequency modes are not restricted to be symmetric 
or antisymmetric any more. Therefore they are coupled 
much more strongly to the bath they lie closer to. 



The disorder changes the normal modes from standing 
waves to localized states. As we are considering disor- 
der in the couplings fi, the high frequency modes are 
affected the most by disorder. We calculate the localiza- 
tion length ^ in units of the lattice constant as the inverse 



3. Temperature profile and heat flux 

Again we transform to real space and calculate the en- 
ergy distribution. There are great differences between the 
single realizations of disorder, so one has to average over 
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(a) 




(b) 





(c) 



FIG. 8: Occupation numbers from 50 realizations of disor- 
der. The lines show Bose-Einstein distributions cor- 
responding to the mean temperature {TaA-Tb) /2, to 
the bath temperatures Ta and Ti, and to the average 
temperature of the normal modes (dashed) . Param- 
eters: Wo = 1, 7 = 2, Ta = 5, Ti, = 2, / = 1 ± 0.2. 
(a) symmetrical chain, Ta = 5, Tt = 2 (b) unsym- 
metrical chain, Ta = 5, Tb — 2 (c) unsymmetrical 
chain, Ta = 0.5, Tt = 0.2 



(Th> 
5 * 



10 20 30 40 50 60 

(a) 



(b) 



FIG. 9: Averaged temperature profile. The error bars show 
the standard deviation of the temperature at each 
lattice site. Parameters: / = 1±0.2, ojq = 1, 7 = 2, 
Ta = 5, Tb — 2. (a) symmetrical, I = 65, 22 realiza- 
tions (b) unsymmetrical, I = 20, 50 realizations 



many realizations in order to obtain comparable results, 
see |Fig. "9} The temperature gradient is enhanced in the 
disordered chain and stays finite, even for long chains. In 
the unsymmetrical case the distribution is wider and the 
average temperature gradient is steeper. In the low tem- 
perature disordered case, |Fig. 10[ the contact resistance 
at the cold bath is large and the temperature gradient 
within the chain remains small. 

<Tr> 
0.75 [ I 



15 



20 



FIG. 10: 



Temperature profile in the unsymmetrical case av- 
eraged over 50 realizations. In the symmetrical 
case the temperature gradient is about half as 
steep. Parameters: / = f ± 0.2, I — 20, ujq = f , 
7 = 2, = 0.5, Tb = 0.2 




10 20 30 40 50 60 70 80 

(a) 



(b) 



FIG. 11: The heat current as a function of the chain length 
I. Each point shows the average over k realizations, 
where 22 < k < 120. The large errorbars show the 
width of the heat flux distribution aj and the small 
errorbars show the error of the average aj/{k — l). 
The solid line is a fit for normal bulk resistance 
proportional to I and the dashed line is a fit for 
a bulk heat resistance proportional to In the 
symmetrical case (a), the chain lengths are large 
enough to indicate that the data fits better to the 
bulk heat resistance proportional to In the 

unsymmetrical case (b) the achieved chain lengths 
are insufficient. Parameters: ljo = 1, 7 = 2, Tq = 
5, Tb = 2, f = 1±0.2. 

The heat flux is reduced in the disordered case, com- 
pare |Fig. 11| and |Fig. 5| (a) . The main difference com- 
pared to the ordered case (see |Fig. 5[ ) is, that the heat 
flux is not independent of the chain length for / > 5 any 
more. In the symmetrical case it was possible to calculate 
the correlations for a chain length up to Z = 75 in a rea- 
sonable time. One can try to determine the asymptotic 
behavior from this data. For heat conduction according 
to Fourier's law we would expect a heat flux proportional 
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to {Rc + -RbuikO""'^ with a contact resistance Rc and the 
resistance of the inner chain i?buik^- For harmonic chains 
the classical asymptotic resistance in proportional to Vl 
0, The dashed lines are fits to this behavior. Due to 
the wide distribution of the heat currents and the finite 
lengths of the considered chains, it is not evident which 
fit describes the asymptotic behavior, but at least in the 
case of symmetric chains our model seems to show 
asymptotic behavior (dashed line). 

C. Entanglement 

Beside the zero point energies and the Bose-Einstein 
statistics of the occupation numbers, entanglement is an- 
other relevant quantum mechanical feature that can be 
observed in our model. 

We use the logarithmic negativity, which is a mea- 
sure for the entanglement between two parts of a system, 
which can be calculated from the correlations between 
the coordinates and momenta of the system [H, [2^ . 

We set up the covariance matrix containing all the cor- 
relations between coordinates and momenta in real space. 
Then the system is divided into two parts A and B and 
the covariance matrix is partially transposed with respect 
to B. The logarithmic negativity is then given by 

7V--^log2(min(l,|7,|)) , (24) 

j 

where the jj are the symplectic eigenvalues of the par- 
tially transposed covariance matrix. We use the nota- 
tion iVfc for the logarithmic negativity of the subsystems 
A = {Xi,...,Xk} a.ndB = {Xk+u ■ ■ ■ , Xi}. Other divi- 
sions, e. g. taking every second oscillator or performing 
the division in normal mode space, are also possible, but 
their physical meaning is not obvious. 

1. Entanglement m dependence of the couplings f and 7 

We start with the ordered case. Figure [T^] shows the 
entanglement for different divisions of a chain of length 
20 as a function of the coupling within the chain /. With 
/ = the oscillators are not coupled at all and there is 
no entanglement. Increasing / favors entanglement. For 
each Nk there is a threshold coupling where entangle- 
ment starts. In iVi and iV;_i one of the subsystems is a 
single oscillator at the end of the chain coupled directly 
to a bath. In this case we observe a lower logarithmic 
negativity and a higher threshold coupling for the onset 
of entanglement. All the other Nk with 1 < k < I — 1 
behave very similarly. 

2. Entanglement as a function of temperature 

We expect the logarithmic negativity to decrease and 
finally vanish with increasing temperature. Additionally 



N 

0.4 



0.3 




FIG. 12: Stationary logarithmic negativity as a function of 
the coupling / inside the ordered chain. Solid line: 
A''i, dashed line: N5. Parameters: I = 20, 7 = 2. 

to the ordered case we want to investigate the entangle- 
ment in the unsymmetrical disordered case. Therefore 
we vary the mean temperature = (Tq -|- Tb)/2 for 
each realization of disorder. Another quantity marking 
the transition from the quantum mechanical regime is the 
heat conductivity (see Sec. lIII A^ which will be observed 
simultaneously. 

Ta and Ti, are chosen as (1 ± e)Tm with e = 0.1. The 
further parameters are chosen 

1^20 r = 10 7 = 2 f = l af^ 0.2 ujQ^l. 

(25) 

In |Fig. 13| (a) the conductivity and entanglements for dif- 
ferent divisions of the chain are plotted for an ordered 
chain. The entanglements with one subsystem consist- 
ing of a single oscillator A^i and A^ig have lower values 
than the others, as observed already in the previous sec- 
tion. All Nk are initially constant and start decreasing 
at T„i ~ 0.2. A^i and A^ig reach zero at w 0.45, the 
others follow at T,,,, « 0.67. In the same temperature 
range the heat conductivity J/ AT rises. 

In the disordered case ( |Fig. 13| (b)) there are only mi- 
nor changes to the entanglement in each realization of 
disorder. The sharp transition to zero is washed out by 
the average. As seen before, the heat conductivity J/ AT 
is lower in the disordered case, but it shows essentially 
the same temperature dependence. 

The plateau we observe at low temperatures results 
from the frequencies of the normal modes starting from 
luq. In the case of a chain with onsite potentials only at 
the ends of the chain, |Fig. 13| (c), the logarithmic neg- 
ativity reaches higher values at low temperatures, but 
their decrease starts already at lower temperatures. 

D. Time evolution 

Finally we shortly want to study some time dependen- 
cies of the system. For simplicity we inspect a short 
ordered chain, consisting of only four oscillators. We 
evaluate the time dependent correlations (like equation 
pop ) and study e. g. the diagonal momentum correla- 
tions, which are proportional to the kinetic energy, see 
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FIG. 13: The heat conductivity J/ AT (dashed hne) and 
some entanglements as a function of temperature. 
Dotted line: A^i; gray line: solid line: N5; 

gray dashed line: Parameters from (I25|l . (a) 

Ordered Chain, e — 0.1 (b) Unsymmetrical dis- 
ordered chain, / = 1 ± 0.2, averaged over 20 re- 
alizations, e = 0.1 (c) Ordered Chain, no onsite 
potential within the chain, e = 0.01 



|Fig. 14| Under some oscillations each lattice site gains 
energy until it reaches the value corresponding to the 
stationary temperature profile. 

Furthermore we generalize equation (|16p for correla- 
tions of two coordinates or momenta at different times. 
In |Fig. 15| results are shown for the time shifted auto- 
correlation functions of the momenta in normal coordi- 
nates in the stationary limit. Each correlation performs 
damped oscillations with its slightly detuned frequency. 
Each normal coordinate loses energy by damping which is 
replaced by incoherent noise. The stronger the influence 
of damping and noise to a normal coordinate, the faster 
it loses memory of its history and the autocorrelation 
function decays. Normal modes with large amplitudes at 
the ends of the chain experience the strongest damping. 




. . . t 

5 10 15 20 



FIG. 14: (Color online) The time evolution of the diagonal 
momentum correlations. Parameters: / = 1, 7 = 
0.5, = 5, Tb = 2 

|ijQ(t),Q<t+r))| 




FIG. 15: (Color online) The time shifted diagonal momen- 
tum correlations in normal coordinates in the limit 
t ^ 00. Parameters like in |Fig. 14\ 

IV. SUMMARY 

We have treated disordered harmonic chains with the 
Quantum Langevin formalism in the strong coupling 
regime. The strong coupling to the heat baths led to 
a renormalization of the normal frequencies. 

In ordered chains the occupation numbers calculated 
with these frequencies approximately follow the Bose- 
Einstein statistics, even in the strong-coupling non- 
equilibrium regime. The heat flux follows from non- 
diagonal correlations between coordinates and momenta. 

Symmetric disorder does not change the occupation 
numbers qualitatively, because the normal coordinates 
are constrained to have the same amplitude at both ends 
of the chain. With breaking the left-right symmetry this 
changes. The localization of most of the modes enhances 
the effect of strongly asymmetric coupling to the heat 
baths. 

In the low temperature regime the energy distribution 
in the chain is dominated by zero-point energies that 
have to be taken into account for constructing the local 
temperature Tr. 

In the limit of ordered chains we have recovered 
the vanishing temperature gradient and the length- 
independent heat flux known from classical models. In 
disordered chains the temperature gradient is finite and 
the heat flux decreases with the length. This decrease 
seems to be slower than l^^, following the classical predic- 
tion for harmonic disordered chains , but the asymp- 
totic behavior could not be identified clearly due to nu- 
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merical restrictions. 

Characteristic quantum mechanical features are the 
freezing of the heat conductivity, that behaves typical for 
optical phonons in our model. In the same temperature 
range where the heat conductivity freezes, entanglement 
appears. 



Appendix A: THE NOISE RESPONSE 
FUNCTIONS IN THE SYMMETRICAL CASE 

1. The structure of the response functions 

We want to calculate only the noise response functions 
and therefore omit the initial conditions in p^ . We use 
the symmetry relations pT|) and get 



jecr(i) 



(Al) 

Dividing this equation by Gu and subtracting the same 
equation with i replaced by n eliminates the 7 and the rj 
terms and yields 



Gin 



Gli 



Using this equation for eliminating all Yj with j ^ i from 
(|XT|) yields 



(A2) 



with 



yields the following form 



{s) = (s + T)Gu n + /^p(^) (^) ■ (A3) 



The numerator contains any (s^ + ^f) from the same 
symmetry, except for the term with its own frequency. It 
is proportional to the amplitude of the mode at the end 
of the chain Gu. 

In the unsymmetrical case this analytic calculation of 
the noise functions is not possible due to mixing of the 
symmetry families. Therefore Cramer's rule has to be 
applied explicitly, which requires a higher numerical ef- 
fort. 

2. The roots of the denominator 

The denominator of the response functions reads 



7r 



We will prove that the real parts of the roots of Dp(i) (s) 
are negative. As s = and s — are no roots, we can 

search for the roots of Dp(^i^{s)/ njeo-(i)(*^ 



which implies 



^2 



(A5) 



Gli 



(s2 + f7f) 



1+E 

ie<T(i) 



2G2 



57(5) 



Extracting the common denominator Z)p(j-)(s), which is 
equal for all response function from the same symmetry, 



Without loss of generality we assume that Im(s) > 0. 
Now we assume Re(s) > 0, which implies arg(s) — 
arg(— 1/s) S (0, 7r/2). On the other hand we find 
arg(s2 + e (0,^) ^ arg(l/(s2 + n^^)) e {-tt,0) 
arg(u(s)) e (— TT, 0). Therefore (|A5p cannot be fulfilled 
and the assumption Re(s) > is disproven. 
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